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Since long back, scientists have been putting enormous effort to understand earthquake dynamics 
-the goal is to develop a successful prediction scheme which can provide reliable alarm that an earth- 
quake is imminent. Model studies sometimes help to understand in some extend the basic dynamics 
of the real systems and therefore is an important part of earthquake research. In this report, we 
review several physics models which capture some essential features of earthquake phenomenon and 
also suggest methods to predict catastrophic events being within the range of model parameters. 



Introduction 



'Earthquake' is a long standing problem to the scientist 
community as it is a subject of many unknowns. Strong 
earthquakes often causes huge devastation in terms of 
lives and properties and therefore this subject always de- 
mands high priority. But so far, little advancement has 
been made to understand the entire earthquake-dynamics 
and the challenge still remains -to develop a successful 
prediction scheme which can provide exact space-time 
information of future earthquakes and their expected 
magnitudes. However, several models and hypothesis 
PL [3- 0, 0, 0, @| have been proposed and some recent 
studies suggest potential methods to predict catastrophic 
events in some earthquake models. In this short report 
we discuss such physics models of earthquake, giving im- 
portance to their capability of prediction. 

Geological facts 



Plate-tectonic theory explains the origin of earthquake 
through a stick-slip dynamics: Earth's solid outer crust 
(about 20 km thick) rests on a tectonic shell which is 
divided into numbers (about 12) of mobile plates, hav- 
ing relative velocities of the order of few centimeters per 
year. This motion of the plates arises due to the power- 
ful convective flow of the earth's mantle at the inner core 
of earth. On the other hand solid-solid frictional force 
arises at the crust-plate boundary and it sticks them to- 
gether. This kind of sticking develops elastic strains and 
the strain energy gradually increases because of the uni- 
form motion of the tectonic plates. Therefore a compe- 
tition comes to play between the sticking frictional force 
and the restoring elastic force (stress). When the accu- 
mulated stress exceeds the frictional force, a slip (earth- 
quake) occurs and it releases the stored elastic energy in 
the form of sound, heat and mechanical vibrations. It 
has been observed that generally a series of small earth- 
quakes appear before (foreshocks) and after (aftershocks) 
a big quake (main shock). 
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Figure. 1: The magnitude distribution of earthquakes in Japan 
(from JUNEC catalog) 0. The straight line with slope -0.9 is 
the best fit to the data points for m > 3 and supports Gutenberg- 
Richter law. 



Physics models 

The overall frequency distribution of earthquakes in- 
cluding foreshocks, mainshocks and aftershocks, seems 
to follow the empirical Gutenberg- Richter law 

In N(m) = constant — bm, 

where N(m) is the number of earthquakes having magni- 
tude (in Richter scale) greater than or equal to m, b is the 
power exponent. The observed value of b ranges between 
0.7 and 1.0 (Fig.l). The amount of energy e released in 
an earthquake is related to the magnitude as 

In e = constant + am. 

Therefore Gutenberg-Richter law can be expressed as an 
alternative form: 

N(e) = e~ a , 

where N(e) is the number of earthquakes releasing energy 
greater than or equal to e and a = b/a. The Gutenberg- 
Richter law is being considered as one of the most fun- 
damental observations by the physicists. 

Several models have been proposed to study the nature 
of the earthquake phenomenon. The main intension is to 
capture the Gutenberg-Richter type power law for the 
frequency distribution of failures (quakes) by modelling 
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different aspects of faults -structure, material properties, 
geometry etc. These models can be classified into four 
groups according to their basic assumptions: (A) Fric- 
tion models, incorporate the stick-slip dynamics through 
the collective motion of an assembly of locally connected 
elements subject to slow driving force (B) Fracture mod- 
els, look at earthquake phenomena as a fracture-failure 
process of deformable materials that break under ex- 
ternal loading through slow build-up of stress (C) Self- 
organised critical (SOC) models, consider earthquake as 
a self-driven slow process and (D) Fractal models, give 
importance to fractal nature of the crust-plate interfaces 
and address the phenomenon as a two-fractal overlap 
problem. 

(A) Friction models 

In 1967 Burridge and Knopoff Q introduced model 
studies in earthquake research. They proposed a spring- 
block model to mimic the typical stick-slip dynamics of 
earthquake phenomena, which has been extended later 
by Carlson and Langer @|. The Burridge-Knopoff type 
model contains a linear array of blocks of mass m coupled 
to each other by identical harmonic springs of strength 
k c and also attached to a fixed surface at the top by a 
different set of identical springs having strength k p . The 
blocks are kept on a horizontal platform (rough surface) 
which moves with a uniform velocity V (Fig. 2). Here 
qualitatively the blocks can be thought of as the points of 
contact between two plates moving at a relative speed V, 
where the spring constants k c and k p represents the linear 
elastic response of the contact region to compression and 
shear respectively. 
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and dynamic phase begins. Due to the presence of fric- 
tional force this dynamic phase is definitely dissipative 
in nature. This dissipation reduces the relative velocity 
of the block and the system again goes back to the static 
phase. The dynamic friction is assumed as a velocity 
weakening function / (shown in the Fig. 2 ). Presence of 
spring force and velocity weakening friction force leads 
the system to a complex dynamical state. At the initial 
stage, individual small slip occurs. But due to constant 
pulling, the springs are gradually stretched and attain 
the limit of total frictional stability of the blocks where 
the collective slips of almost all the blocks occur. Clearly 
this is a big event and a large amount of elastic energy 
is released here. The equation of motion of the j th block 
of the system is 



k c [x 



l) 



f(±j - V) 



where dots denote differentiation with respect to t, m is 
the mass of a block, k p and k c are spring constants of the 
connecting springs and / represents the nonlinear veloc- 
ity weakening friction force. The position coordinate x 
is measured along the chain length and the velocity V of 
the platform is in the increasing x direction. The total 
energy of the system at time t can be calculated by solv- 
ing the above equation. A sudden drops in the energy 
of the system is identified as the released energy during 
a slip event and the distribution of such energies follows 
power-law: 



N(e) 



1, 



where N(e) is the number of slips releasing energy greater 
than or equal to e. 

Prediction possibility of major events (slips) 

Recently Dey et al. [§] developed a method to predict 
the major slip event in Burridge-Knopoff type spring- 
block model. Introducing an additional dissipative force 
in the spring-block arrangement they identified the dis- 
sipative functional R(t) as energy bursts similar to the 
acoustic emission signals ficl fill ] observed in experi- 
ments. The distribution of R(t) shows power laws if one 
records all slip events including the major slips. 



Figure 2: The Burridge-Knopoff model (above) and two differ- 
ent forms of velocity weakening friction force (below) [jj. 

Starting from unstrained condition, when the system is 
pulled slowly (at constant rate), the blocks initially re- 
main stuck to the surface due to friction between sur- 
faces. A slip of a block occurs when the correspond- 
ing spring force overcomes the threshold value (maxi- 
mum static frictional force between that block and the 
rough surface). The dynamics of each block is basically 
the resultant of two phases: a static phase and a dy- 
namic phase. During static phase the block remains stuck 
on the rough surface and the elastic strain continuously 
grows up. The static phase comes to a sudden end when 
driving force on that block attains the threshold value 




Figure 3: The dissipative functional R(t) vs. time (t) in 
Burridge-Knopoff model [3|. 

The plot of R(t) vs. t shows a gradual increase in 
activity (Fig. 3) prior to the occurrence of a major slip. 
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But as R(t) is noisy it can not help much to predict the 
major slip event, rather the cumulative energy dissipated 
Eaeit) ~ Jo R(t)dt grows in steps and it seems to di- 
verge as a major slip event is approached (Fig. 4). From 
such divergence one can predict the occurrence time (t c ) 
of a major slip through proper extrapolation. 
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Figure 4: Cumulative energy E ae versus time (t) plot. Inset 
shows the time series of R(t) including two major slip events [jj. 

(B) Fracture models 



Earthquake can be considered as a fracture-failure phe- 
nomenon through slow build up of stress at the plate- 
crust boundary with the movement of the plate as exter- 
nal driving force. The deformation properties of the ma- 
terials sitting at the boundary play crucial role on spatial 
redistribution of the stress around a broken region and 
this actually decides in which direction the crack front 
should propagate. Fiber bundle model and random fuse 
model address such scenario in material breakdown and 
sometimes they are treated as models of earthquake since 
they produce time series of avalanches which follow power 
law distribution similar to Gutenberg-Richter law. 



Fiber bundle model 

Fiber bundle (RFB) model consists of many (No) fibers 
connected in parallel to each other and clamped at their 
two ends and having randomly distributed strengths. 
The model exhibits a typical relaxational dynamics when 
external load is applied uniformly at the bottom end 
(F i g.5) . In the global load-sharing approximation fl2l . 
LL3J, LLJ], surviving fibers share equally the external load. 
Initially, after the load F is applied on the bundle, fibers 
having strength less than the applied stress a = F/Nq 
fail immediately. After this, the total load on the bun- 
dle redistributes globally as the stress is transferred from 
broken fibers to the remaining unbroken ones. This re- 
distribution causes secondary failures which in general 
causes further failures and produces an "avalanche" which 
denotes simultaneous failure of several elements. With 
steady increase of external load, avalanches of different 
size appear before the global breakdown where the bun- 
dle collapses. The scaling properties of such mean-field 



dynamics and the avalanche statistics are expected to be 
extremely useful in analysing fracture and breakdown in 
real materials, including earthquakes [H, EE G3, E| ■ 



Figure 5: The fiber bundle model . 

When external load is applied, the surviving fraction 
of total fibers follows a simple recursion relation 

U t+1 (a) = 1-P(a/U t ), 

where Ut = Nt/No, a is the external stress and P is 
cumulative probability function. The recursion relation 
has the form of an iterative map Ut+i = Y(Ut) and finally 
the dynamics stops at a fixed point where Ut+i = Ut. 

If external load is increased in steps by equal amount 
AF, then the entire failure process can be formulated 
through the recursive dynamics mentioned above and 
the fixed point solution gives the value of critical stress 
a c above which the bundle collapses. It can be shown 
that the bundle undergoes a phase transition from par- 
tially broken state to completely broken state. The order 
parameter (O), susceptibility (x) and relaxation time (r) 
follow robust power laws with universal exponent values 

In case of quasi-static load increment, only the weakest 
fiber (among the intact fibers) fails after loading and then 
the bundle undergoes load redistribution till a fixed point 
is reached. The fluctuations in strength distributions pro- 
duces different size of avalanches during the entire failure 
process and the avalanche distribution follows power law 
(Fig. 6) with exponent —5/2. Hemmer and Hansen 
has analytical proved that this exponent is universal un- 
der mild restriction on strength distributions. 
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Figure 6: Avalanche time series in a fiber bundle model (left) 
and the corresponding avalanche distribution (right). 
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Random fuse model 




Figure 7: The random fuse model. 

The fuse model [U, H2, H3] consists of a lattice in which 
each bond is a fuse, i.e., an ohmic resistor as long as the 
electric current it carries is below a threshold value. If the 
threshold is passed, the fuse burns out irreversibly. The 
threshold t of each bond is drawn from an uncorrelated 
distribution p(t). The lattice is placed at 45° with re- 
gards to the electrical bus bars (Fig. 7) and an increasing 
current is passed through it. Numerically, the Kirchhoff 
equations are solved at each node. 

When a bond breaks, current value on the neighbor- 
ing bonds increases and sometimes it triggers secondary 
failures. Finally the system comes to a state where no 
current passes through the lattice; that means there is 
a crack which separates the lattice in two pieces. With 
gradual increase of current / voltage a series of intermedi- 
ate avalanches appear before the final breakdown. The 
distribution of such avalanches follows power law with 
exponent close to 3 (Fig. 8). 
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Figure 9: Crossover signature in avalanche power law. In fiber 
bundle model: xq is the starting position of recording avalanches 
and x c is the global failure point; exponent value changes from 5/2 
to 3/2 (left) and in random fuse model: exponent value changes 
from 3 to 2 when the system comes closer to breakdown point 
(right). 

Recently Kawamura [27] observed similar crossover 
behavior for the local magnitude distribution of earth- 
quakes in Japan (Fig. 10). This observation has strength- 
ened the possibility of using crossover signal as a tool of 
predicting catastrophic events. 
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Figure 8: Avalanche time series in a fuse model (left) and the 
corresponding avalanche distribution (right). 



Crossover behavior: signature of imminent breakdown 

A robust crossover behavior has been observed [13, US 
|26| recently in the two very different models described 
above where the system gradually approaches the global 
failure through several intermediate failure events. If in- 
termediate avalanches are recorded, avalanche distribu- 
tion follows a power law with an exponent that crosses 
over from one value to a very different value when the 
system is close to the global failure or breakdown point 
(Fig. 9). Therefore, this crossover is a signature of immi- 
nent breakdown. 



Figure 10: Crossover si gna ture in the magnitude distribution of 
earthquakes within Japan [27l |. 

(C) Self-organised critical models 

In various thermodynamics systems, there is a "criti- 
cal" point where the systems become totally correlated 
and show scale free (power law) behaviour. Apart from 
the critical point the average microscopic quantities of 
the systems follow scaling behaviour. Generally such crit- 
ical states are achieved through the fine-tuning of phys- 
ical parameters, such as temperature, pressure etc. and 
the power law behavior is considered to be a signature of 
the "critical" state of the system. However, it has been 
observed that some complex systems evolve collectively 
to such critical state only through mutual interactions 
and show power law behavior there. These systems do 
not need any fine tuning of physical parameter and there- 
fore are considered as "self-organised critical" (SOC) sys- 
tems. The term SOC was first introduced by Bak, Tang 
and Wiesenfeld in 1987. 

The magnitude distribution of earthquake shows 
power-law (Gutenberg- Richter law), therefore it is tempt- 
ing to assume that earthquake happens through a self- 
organised dynamics: the build up of stress due to tectonic 
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motion is a self organised slow process; gradually the crit- 
ical state is achieved where the stress releases in bursts 
of various sizes. Extensive research have been going on 
to establish relation between earthquake and SOC sys- 
tems, for which several models have been proposed. So 
far, sandpile models are the best example of SOC system. 



ultimately settles at a critical state having height h c and 
exhibits scale free behavior in terms of avalanche and 
life time distributions. But the power law exponents are 
different compared to those in BTW model and there- 
fore Manna model belongs to a different universality class 




BTW sandpile model 



The first attempt to study SOC through model systems 
was made by Bak, Tang and Wiesenfeld This is 

a model of sandpile whose natural dynamics drives it 
towards the critical state. The model can be described 
on a two dimensional square lattice. At each lattice site 
there is an integer variable h%,j which represents 
the height of the sand column at that site. A unit of 
height (one sand grain) is added at a randomly chosen 
site at each time step and the system evolves in discrete 
time. The dynamics starts as soon as any site (i, j) has 
got a height equal to the threshold value (h t h= 4): that 
site topples, i.e., hij becomes zero there, and the heights 
of the four neighbouring sites increase by one unit 



Critical state 



hi 



4, fk±x t j — > hi±i,j + 1, hij±i — > hij±i + 1. 
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Figure 12: Avalanche time series in a BTW model (left) and in 
a Manna model (right). 
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Figure 11: BTW model on square lattice [29l |. 

The process continues till all sites become stable 
(Fig. 11). In case of toppling at the boundary of the 
lattice (4 nearest neighbours are not available), grains 
falling outside the lattice are removed and considered 
to be absorbed/collected at the boundary. Gradually 
the average height h av attains a critical value h c , be- 
yond which it does not increase at all - on an average 
the additional sand grains are flown away from the lat- 
tice . Total number of toppling between two successive 
stable states, determines the size of an avalanche and at 
the critical state avalanche size distribution follows power 
laws: n s ~ s~ r , where n s denotes the density of s size 
avalanches. The exponent F has the value T ~ 1.15±0.10 
in 2D 1291. 



Manna model 

Manna proposed the stochastic sand-pile model [2§| 
by introducing randomness in the dynamics of sand-pile 
growth. Here, the critical height is 2. Therefore at each 
toppling, the two rejected grains choose their host among 
the four available neighbours randomly with equal prob- 
ability. After constant adding of sand grains, the system 



Precursory activities 

On constant adding of grains, sandpile models gradu- 
ally attains critical state from sub-critical states. Now 
the question is: can one predict the critical state from 
the sub-critical response of the pile? A pulse perturba- 
tion method [U, [32] gives the answer: 

At an average height h av , a fixed number of height 
units hp (pulse of sand grains) is added at any central 
point of the system. Just after this addition, the local 
dynamics starts and it takes a finite time or iterations 
to return back to the stable state after several toppling 
events. One can measure the response parameters: A — > 
number of toppling r — » number of iteration and £ — > 
correlation length which is the distance of the furthest 
toppled site from the site where h p has been dropped. 

To ensure toppling at the target site, the pulse hight 
has been chosen as h p — 4 for BTW model and h „ = 
2 for Manna model. Simulation studies conclude [32] 
that all the response parameter follow power law as h c 
is approached: A oc (h c — h av )~ x , r oc (h c — h av )~^, 
i oc (h c - h a v)~ v ] A £3 2.0, /i = 1.2 and v S 1.0. Now 
if A" 1 / A , t- 1 ^ and f -1 /" are plotted against h av , all 
the curve follow straight line and they should touch the 
x axis at h av = h c . A proper extrapolation estimates 
the critical hight h c = 2.13 ± .01 for BTW model and 
h c = 0.72 ± .01 for Manna model, which agree well with 
direct estimates I28B. 




Figure 13: Sub-critical response: Correlation length (£) versus 
average hight (h av ) in BTW model (left) and in Manna model 
(right). Inset shows the plot of inverse £ versus h av and predicts 
the critical hight through extrapolation. 

Therefore, although BTW and Manna models belong 
to different universality classes with respect to their prop- 
erties at the critical state, both the models show similar 
sub-critical response or precursors. A proper extrapola- 
tion method can estimate the respective critical heights 
of the models quite accurately. 

(D) Fractal models 

The surfaces of earth's crust and tectonic plate at the 
fault zone are not compact rather fractal in nature. In 
fact, these surfaces are the results of the large scale frac- 
ture separating the crust from the moving tectonic plate. 
It has been observed that these surfaces are self sim- 
ilar fractals [33| having the self-affine scaling property 
h(\x,\y) ~ X^h(x,y), where h(x) denotes the height of 
the crack surface at the point x and ( is the roughness 
exponent. It has been claimed recently that since the 
fractured surfaces have got well-characterized self-affine 
properties, the distribution of the elastic energies released 
during the slips (earthquake events) between two rough 
surfaces (crust and plate) may follow the overlap distri- 
bution of two fractal surfaces (34I l35| 




Figure 14: A typical self-affine fracture surface in (2 + 1) dimen- 
sion with roughness exponent £ = 0.8. 



Self-affine asperity model 

V. De Rubies et. al. (34| has proposed a new model for 
earthquakes where the scale invariance of the Gutenberg- 
Richter law is claimed to come from the fractal geome- 
try of the fault surfaces. In this model the sliding fault 
surfaces have been represented by fractional Brownian 
surfaces, whose height scales as \h(x + r) — h{x)\ ~ r 1 *. 
The roughness of the surfaces are determined from the 
value of the exponent C, which lies between and 1. Two 
surfaces are simulated by two statistically self-affine pro- 
files, say, h\{x) and fi2(x), one drifting over other with a 
constant speed v such that hi(x, t) = hi(x — vt). An in- 
teraction between two profiles represents a single seismic 
event and the energy released is assumed to be propor- 
tional to the breaking area of the asperities. At the con- 
tact point of the surfaces the 'surface roughness' prevents 
slipping and the stress is accumulated there. When the 
stress exceeds a certain threshold value, breaking (earth- 
quake) occurs. This model produces Gutenberg-Richter 
type power law: P(E) ~ E^ ' 1 , where P{E)dE is the 
probability that an earthquake releases energy between 
E and E + dE and the exponent j3 is directly related 
to the roughness exponent ( &s : (3 = 1 — (/{d — 1) = 
(Dj — \)/{d — 1), where Df and d are respectively the 
fractal dimension and the embedding dimension of the 
surfaces. 



Chakrabarti-Stinchcombe model 

This is an analytical model, proposed by Chakrabarti and 
Stinchcombe [35||, which incorporates the self-similar na- 
ture of both the crust and the tectonic plate. They used 
self-similar fractals to represent fault surfaces (Fig. 15). 




u u n n 



Figure 15: Schematic representation of the rough surfaces of 
earth's crust and moving tectonic plate [35| . 

The total contact area between the surfaces is assumed 
to be proportional to the elastic strain energy that can 
be grown during the sticking period, as the solid-solid 
friction force arises from the elastic strain at the con- 
tacts between the asperities. This energy is considered 
to be released as one surface slips over the other and 
sticks again to the next contact between the rough sur- 
faces. Chakrabarti and Stichcombe have shown analyti- 
cally through renormalization group calculations that for 
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regular fractals (Cantor sets and carpets) the contact 
area follows power law distribution: 



P(s) 



;7 = i. 



which is comparable to that of Gutenberg-Richter law. 



Numerical verification 

The claim of Chakrabarti-Stichcombe model has been 
verified by extensive numerical simulations [3(| taking 
different type of synthetic fractals: regular or non- 
random Cantor sets, random Cantor sets (in one dimen- 
sion), regular and random gaskets on square lattice and 
percolating clusters embedded in two dimensions. 




Figure 17: The plot of P'(m') against m' for Cantor sets with 
df = In 2/ In 3 at finite generations: n = 10 (square), n = 11 (plus), 
n = 12 (cross) and n = 13 (star) and the dotted lines indicate the 
best fit curves of the form a(x — b)~ d ; where d = 1. Inset shows 
P(m) vs. m plots. 




Figure 16: (a) A regular Cantor set of dimension In 2/ In 3; only 
three finite generations are shown, (b) The overlap of two identical 
(regular) Cantor sets, at n = 3, when one slips over other; the 
overlap sets are indicated within the vertical lines, where periodic 
boundary condition has been used. 



The contact area distributions P(m, L) seem to follow 
a universal scaling: 

P(m, L) ~ L a P'(m')-m' = mL a , 

where L denotes the size of the fractal and a = 2(d — 
df); df being the mass dimension of the fractal and d is 
the embedding dimension. Also the overlap distribution 
P(m), and hence the scaled distribution P'(m'), decay 
with m or m' following a power law (Fig. 17) for both 
regular and random Cantor sets and gaskets: 

P(m) = m- f3 -j3 = d. 

A very recent report [13] analytically explains the ori- 
gin of such asymptotic power laws in case of Cantor set 
overlap. 



Prediction possibility of large events 

If one cantor set moves uniformly over other, the over- 
lap between the two fractals change quasi-randomly with 
time and produces a time series of overlaps m(t). Such a 
time series is shown in Fig. 18, for Cantor sets of dimen- 
sions In 2/ In 3. While most of the overlaps are small in 
magnitude, some are really big, where as the cumulative 
overlap size Q{t) = J mdt 'on average' grows linearly 
with time. Is it possible to predict a large future over- 
lap analysing the time series data? A recent study [H| 
suggests a method: 



, — .150 - 




Time (t) 



Figure 18: The time (t) series data of overlap size (m) for 
regular Cantor sets: of dimension In 2/ In 3, at 8th generation. 

One can identify the 'large events' occurring at time ti 
in the m(t) series, where m(ti) > M, a pre-assigned num- 
ber, then calculate the cumulative overlap size Q(t) = 
J t * +1 mdt, where the successive large events occur at 
times ti and t i+ i. Obviously Q{t) is reset to value 
after every large event. The behavior of Qi with time is 
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shown in Fig. 19 for regular cantor sets . It appears that 
there are discrete values up to which Qi grows with time. 




Time (tj Time (t) 



Figure 19: The cumulative overlap size variation with time 
(for regular Cantor sets of dimension In 2/ In 3, at 8th genera- 
tion), where the cumulative overlap has been reset to value 
after every big event (of overlap size > M where M = 128 
and 32 respectively). 



Therefore, if one fixes a magnitude M of the overlap 
sizes m, so that overlaps with m> M are called 'events' 
(or earthquake), then the cumulative overlap Qi grows 
linearly with time up to some discrete levels Qi = IQo, 
where Qq is the minimal overlap size, dependent on M 
and I is an integer. This information certainly does not 
help to predict a future large event accurately, but it gives 
some hints by identifying discrete levels of Qi where a 
large overlap is likely to happen. 

Discussions and concluding remarks 

The Gutenberg-Richter law is a well established law in 
earthquake research. It should be mentioned that Guten- 
berg and Richter obtained this law from the statistics of 
earthquake events observed throughout the world. Also, 
earthquakes within a tectonically active region (Japan, 
California etc) follow similar power law. It is still a con- 
troversial issue whether the exponent of the Gutenberg- 
Richter power law is an universal constant or it varies 
in a narrow range (0.8 to 1.2) depending upon the na- 
ture of the fault zone. As the motion of the tectonic 
plates is surely an observed fact, the stick-slip process 
should be a major ingredient of earthquake models. Al- 
though Burridge-Knopoff type spring-block model suc- 
cessfully captures the stick-slip dynamics and reproduces 
Gutenberg-Richter type power law in the size distribu- 
tion of events, it is far from the accurate representation of 
earthquake dynamics -as it does not contain a mechanism 
for aftershocks which are observed facts. But this model 
has some important features: The dynamics is inherently 



chaotic -therefore technically unpredictable which agrees 
well with the occurrence of earthquakes. However, the 
method proposed by Dey et. al [3] to predict a major 
slip event by monitoring the cumulative energy function, 
may open up a wide scope of future research in this field. 
On the other hand, fiber bundle model and fuse model 
have been developed basically to study breakdown phe- 
nomena in composite materials. As earthquake is a ma- 
jor breakdown phenomenon, these models can be used 
as earthquake models due to their inherent mean-field 
nature. Avalanche distributions show power laws in both 
the models and a crossover in exponent value appears [25| 
near breakdown point, which can be treated as a crite- 
rion for imminent breakdown. SOC models assume a self- 
driven slow dynamics and reproduces Gutenberg-Richter 
law at the critical point. Although the critical state can 
be predicted accurately from the sub-critical response of 
the systems, the behavior remains unpredictable at the 
critical state. Fractal overlap models are different (from 
the other three types) modelling approaches in the sense 
that they focus on the fractal nature of the fault inter- 
faces and not on the dynamics. The contact area distri- 
bution follows asymptotic power law and this suggests a 
possibility that fractal geometry of the faults might be 
the true origin of Gutenberg-Richter law. 

There are several difficulties in studying earthquake 
phenomenon, as it is a "N body" complex problem. While 
exact solution of a 3-body problem needs rigorous math- 
ematical calculations and it does not work for a mutu- 
ally interacting system having more than 3 bodies, nat- 
urally, theoretical physics can not help to formulate the 
entire earthquake dynamics. Again, as the dynamics of 
earthquake happens at a depth more than 20 km from 
earth surface -experimental observation of such dynam- 
ics is almost impossible. Moreover, earthquake dynamics 
involves crust and plates which are highly heterogeneous 
at many scales from the atomic scale to the scale of tec- 
tonic plates, with the presence of dislocation, impurities, 
grains, water etc and the scenario becomes very much 
complicated. In this situation, model studies are very im- 
portant in earthquake research in the sense that they can 
produce synthetic earthquake events, allow us to monitor 
the dynamics and analyse the event statistics to compare 
with the real earthquake data. Moreover, such studies 
suggest potential methods to predict a major event. It 
will be a real breakthrough if any of such methods can 
help a little to predict a future earthquake. 
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